Assignment 7

Lorenzo Biasi, Julius Vernie



In [1]:

    
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from numpy.linalg import inv
%matplotlib inline

We load the variables and initilize the parameters we need



In [2]:

    
data = loadmat('data_files/Tut7_file1.mat')
locals().update(data)
data.keys()









    Out[2]:





dict_keys(['__header__', '__version__', '__globals__', 'A', 'B', 'C', 'Gamma', 'L0', 'Sigma', 'mu0', 'u', 'x', 'z'])



In [3]:

    
p, T = z.shape



In [4]:

    
mu = np.zeros(z.shape)
K = np.zeros((4, 4, T))
V = np.zeros((4, 4, T))
L = np.zeros((4, 4, T))

K[...,0] = L0.dot(B.T.dot(inv(B.dot(L0.dot(B.T)) + Gamma)))
mu[..., [0]] = A.dot(mu0) + K[..., 0].dot(x[:, [0]] - B.dot(A.dot(mu0))) + C.dot(u[..., [0]])
V[..., 0] = (np.eye(4) - K[..., 0].dot(B)).dot(L0)
L[..., 0] = A.dot(V[..., 0].dot(A.T)) + Sigma

We run the filter



In [5]:

    
for t in range(1, T):
    K[...,t] = L[..., t - 1].dot(B.T.dot(inv(B.dot(L[..., t - 1].dot(B.T)) + Gamma)))
    mu[..., [t]] = A.dot(mu[..., [t-1]]) + K[..., t].dot(x[:, [t]] - B.dot(A.dot(mu[..., [t-1]]))) + C.dot(u[..., [t]])
    V[..., t] = (np.eye(4) - K[..., t].dot(B)).dot(L[..., t-1])
    L[..., t] = A.dot(V[..., t].dot(A.T)) + Sigma

We can see a slight offset, we would expect that to be solved with the smoother step



In [6]:

    
plt.plot(mu.T)
plt.plot(z.T, color='red')









    Out[6]:





[<matplotlib.lines.Line2D at 0x7f53ad3f69e8>,
 <matplotlib.lines.Line2D at 0x7f53ab5f9f60>,
 <matplotlib.lines.Line2D at 0x7f53ab603160>,
 <matplotlib.lines.Line2D at 0x7f53ab603320>]



In [7]:

    
V_tilde = np.zeros(V.shape)
mu_tilde = np.zeros(mu.shape)
V_tilde[..., -1] = V[..., -1]
mu_tilde[..., [-1]] = mu[..., [-1]]

for t in range(T - 2, -1, -1):
    #print(t)
    W = V[..., t].dot(A.T.dot(inv(L[..., t])))
    V_tilde[..., t] = V[..., t] + W.dot(V_tilde[..., t+1] - L[..., t]).dot(W.T)
    mu_tilde[..., [t]] = mu[..., [t]] + W.dot(mu_tilde[..., [t+1]] - A.dot(mu[..., [t]]))

we can see that the offset is still present and slightly worse



In [8]:

    
plt.plot(mu_tilde.T)
plt.plot(z.T, color='red')









    Out[8]:





[<matplotlib.lines.Line2D at 0x7f53ab56e908>,
 <matplotlib.lines.Line2D at 0x7f53ab4e5518>,
 <matplotlib.lines.Line2D at 0x7f53ab4e56d8>,
 <matplotlib.lines.Line2D at 0x7f53ab4e5898>]

The predition is clearly following the data more or less correctly, but there is a probelm with the offset, that makes $\tilde{\mu}$ worse than our $\mu$. This should not happen, we would rather expect the opposite.



In [9]:

    
print ('Non smoothed result:', np.sum((mu - z).T ** 2))
print('Smoothed result:', np.sum((mu_tilde - z).T ** 2))

print('Ratio, \n', np.sum((mu_tilde - z).T ** 2) / np.sum((mu - z).T ** 2))









    



Non smoothed result: 9175009.49339
Smoothed result: 9480270.18934
Ratio, 
 1.03327088611

After checking the algorithm many times i decided to look at our x to see if there was anything strange. And if you look closely at the first time steps there is some oddity.



In [10]:

    
plt.plot(x.T)









    Out[10]:





[<matplotlib.lines.Line2D at 0x7f53ab4a1400>,
 <matplotlib.lines.Line2D at 0x7f53ab4a15c0>,
 <matplotlib.lines.Line2D at 0x7f53ab4a17b8>,
 <matplotlib.lines.Line2D at 0x7f53ab4a19b0>]

If someone looks at how x varies at the first time step you will see that it is almost constant and than it starts changing. this could explain the offset in our predictions.



In [11]:

    
#plt.plot(x.T[:4, :])
plt.plot(np.diff(x[..., :10]).T)
np.diff(x[..., :4])









    Out[11]:





array([[   0.7255109 ,  157.7771018 ,  161.32022023],
       [   1.56485095,   52.69063853,   55.13791943],
       [   0.82161329,  -11.52514839,  -13.04291285],
       [   1.33883859,  332.90369003,  340.26390644]])

To test my hunch I decided to remove one time step from the data, to make sure that $x_1$ was not used in the prediction.



In [12]:

    
T = 99
z = z[:, :-1]
mu = np.zeros(z.shape)
K = np.zeros((4, 4, T))
V = np.zeros((4, 4, T))
L = np.zeros((4, 4, T))

K[...,0] = L0.dot(B.T.dot(inv(B.dot(L0.dot(B.T)) + Gamma)))
mu[..., [0]] = mu0
V[..., 0] = 0
L[..., 0] = L0



In [13]:

    
for t in range(1, T):
    #print(t)
    K[...,t] = L[..., t - 1].dot(B.T.dot(inv(B.dot(L[..., t - 1].dot(B.T)) + Gamma)))
    mu[..., [t]] = A.dot(mu[..., [t-1]]) + K[..., t].dot(x[:, [t + 1]] - B.dot(A.dot(mu[..., [t-1]]))) + C.dot(u[..., [t]])
    V[..., t] = (np.eye(4) - K[..., t].dot(B)).dot(L[..., t-1])
    L[..., t] = A.dot(V[..., t].dot(A.T)) + Sigma



In [14]:

    
plt.plot(mu.T)
plt.plot(z.T, color='red')
np.sum((mu - z)**2)
A.dot(mu[..., [t-1]]) + K[..., t].dot(x[:, [t + 1]] - B.dot(A.dot(mu[..., [t-1]]))) + C.dot(u[..., [t]])









    Out[14]:





array([[  -7.65758593],
       [-230.44138841],
       [-592.64195212],
       [-523.77934228]])



In [15]:

    
V_tilde = np.zeros(V.shape)
mu_tilde = np.zeros(mu.shape)
V_tilde[..., -1] = V[..., -1]
mu_tilde[..., [-1]] = mu[..., [-1]]

for t in range(T - 2, -1, -1):
    W = V[..., t].dot(A.T.dot(inv(L[..., t])))
    V_tilde[..., t] = V[..., t] + W.dot(V_tilde[..., t+1] - L[..., t]).dot(W.T)
    mu_tilde[..., [t]] = mu[..., [t]] + W.dot(mu_tilde[..., [t+1]] - A.dot(mu[..., [t]]))



In [16]:

    
plt.plot(mu_tilde.T)
plt.plot(z.T)









    Out[16]:





[<matplotlib.lines.Line2D at 0x7f53ab2a57f0>,
 <matplotlib.lines.Line2D at 0x7f53ab3a8400>,
 <matplotlib.lines.Line2D at 0x7f53ab3a87f0>,
 <matplotlib.lines.Line2D at 0x7f53ab3a8e48>]



In [17]:

    
print ('Non smoothed result:', np.sum((mu - z).T ** 2))
print('Smoothed result:', np.sum((mu_tilde - z).T ** 2))

print('Ratio, \n', np.sum((mu_tilde - z).T ** 2) / np.sum((mu - z).T ** 2))









    



Non smoothed result: 28182.6775642
Smoothed result: 21398.1349435
Ratio, 
 0.759265506082

We can see how by making this slight change to the data the predition are more accurate, and do not have strange sistematical errors.



In [ ]: